RKKY interaction in heavily vacant graphene 
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Dirac electrons in clean graphene can mediate the interactions between two localized magnetic moments. 
The functional form of the RKKY interaction in pristine graphene is specified by two main features: (i) an 
atomic scale oscillatory part determined by a wave vector Q connecting the two valleys. Furthermore with 
doping another longer range oscillation appears which arise from the existence of an extended Fermi surface 
characterized by a single momentum scale kp. (ii) R^ decay in large distances where the exponent a = —3 is 
a distinct feature of undoped Dirac sea (with a linear dispersion relation) in two dimensions. In this work, we 
investigate the effect of a few percent vacancies on the above properties. Depending on the doping level, if the 
chemical potential lies on the linear part of the density of states, the exponent a remains close to —3. Otherwise 
a reduces towards more negative values which means that the combined effect of vacancies and the randomness 
in their positions makes it harder for the carriers of the medium to mediate the magnetic interaction. Addition 
of a few percent of vacancies diminishes the atomic scale oscillations of the RKKY interaction signaling the 
destruction of two-valley structure of the parent graphene material. Surprisingly by allowing the chemical 
potential to vary, we find that the longer-range oscillations expected to arise from the existence of a kp scale in 
the vacant graphene are absent. This may indicate possible non-Fermi liquid behavior by "alloying" graphene 
with vacancies. The complete absence of oscillations in heavily vacant graphene can be considered an advantage 
for applications as a uniform sign of the exchange interaction is desirable for magnetic ordering. 

PACS numbers: 71.23.-k, 73.22.Pr, 71.55.-i, 71.10.Hf 
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INTRODUCTION 

Two-dimensional nature of graphene, along with the 
Dirac nature of charge carriers dll -as contrasted to the 
Schroedinger nature of the carriers in ordinary conductors - 
makes graphene a spectacular platform for condensed matter 
realization of many exciting ideas of low-dimensionality and 
relativistic physics ||2i |3|]. Among the properties of interest 
in materials is the nature of effective interaction between ex- 
ternal agents mediated by the carriers of the host materials. 
An important example of this would be the Ruderman-Kittel- 
Kasuya-Yosida (RKKY) interaction which is the exchange in- 
teraction between two magnetic impurities mediated via the 
propagators of the host. The RKKY interaction in graphene 
was studied first by Saremi [4] who found a generic "atomic- 
scale" oscillatory behavior where the sign of the interaction 
alternates between the two sub-lattices of the bipartite honey- 
comb structure. The wave vector characterizing these oscilla- 
tions is Q = K^ — K~ where K^ denote the momenta asso- 
ciated with two Dirac cones. The oscillations decay as R~^ in 
long distances, where the exponent a = — 3 comes from the 
linear Dirac cone dispersion in two dimensions. The investi- 
gation of Saremi was extended by Sherafati and Satpathy in 
two directions. First for the undoped graphene they extended 
the results to a model of full 7r-bands 1151] . Secondly for the 
low-energy Dirac cone model, they considered the effect of 
non-zero doping ||6|], and as in ordinary two-dimensional con- 
ductors additional /c^ -related oscillations with a R~'^ decay at 
large distances appeared. Such long-wave-length oscillations 
are due to an underlying sharp Fermi surface and disappear 
when the limit of undoped graphene (kp -^ 0) is approached 
in pristine graphene. 

The above line of research was also extended to the case of 



disordered graphene by Lee and coworkers who considered 
the effect of diagonal and hopping disorder on the RKKY 
interaction, in undoped [TJ and doped graphene [S!]. They 
found that for strong enough on-site disorders where Ander- 
son localization takes place, the RKKY interaction exponen- 
tially decays with distance which is consistent with the picture 
based on localized spectrum Jsl] . For the doped graphene and 
weak disorder they found that the Friedel oscillations due to 
the Fermi surface co-exists with the atomic-scale oscillations 
due to interference between the two Dirac cones. In this work 
guided by the experimentally conceivable conditions, instead 
of a generic Anderson model, we consider a specific model 
for vacancies and investigate the RKKY interaction in pres- 
ence of randomly distributed vacancies. Our motivation is 
that various forms of defects are ubiquitous in the surface of 
graphene lIlQII . The presence of defects modifies electronic, 
magnetic and mechanical properties of materials. For exam- 
ple vacancies generated by ion irradiation of graphene give 
rise to magnetic moments Illllll2|] and the most probable form 
of defects produced by irradiation are single vacancies [Il3ll . 
Moreover coherence exchange of spin between such magnetic 
moments and the spectrum of fermions of the ho st graph ene 
leads to unconventional Kondo effect in graphene |ll4i ll5|] . 

In this work we find that the randomness caused by random 
vacancies has a different physics from the Anderson model 
studied in Ref . ||7|] : We find that propagation of electron waves 
in defective environment reduces the exponent a in R'^ dis- 
tance dependence of the RKKY interaction from a = —3 to 
below —4 away from half-filling. The exponent a will not 
differ from —3 as long as the chemical potential is around 
energies where the density of states is linear in energy. We 
find that a is controlled by the concentration of vacancies and 
the chemical potential, but not much dependent on the energy 




FIG. 1. (Color online) Schematic picture of the electronic state 
around the defect in graphene. The symbol at site 1 shows the ac- 
tive sp^ orbital and the circles at sites 2 and 3 the tt orbitals. Two 
of the three sp^ orbitals (i = 2, 3) form a covalent bond. V is the 
amplitude of the hybridization between the active sp^ orbital and the 
two TT orbitals. 



level of the localized orbitals bound to vacancy location. At 
higher concentrations of vacancies, they are expected to form 
their own impurity band, which is expected to cause a metal- 
lic type transport in neutral vacancies when the concentration 
goes beyond a few percent |tl6i \V^ . In such case we expect a 
Friedel oscillations to appear as a signature of metallic state. 
But surprisingly we find that in addition to atomic-scale os- 
cillations, the Friedel oscillations of the Fermi surface also 
disappear. 



MODEL AND METHOD 

Kanao and coworkers proposed a simple and conceivable 
model for a localized orbital and its hybridization with the ad- 
jacent Pz orbitals needed to understand the physics of local 
magnetic moments and their interaction with the Dirac elec- 
trons of graphene [18]. They considered a Jahn-Teller dis- 
torted geometry of a single atom point defect shown in Fig. [T] 
where the sp^ orbital of one of the atoms, e.g. 1 is hybridized 
with Pz orbitals of the two remaining atoms 2,3 surround- 
ing the vacancy. The effective Hamiltonian for the defective 
graphene becomes llisll : 



H = Hg-\- i/def + ^hyb , 



(1) 



where Hg is the tight-binding Hamiltonian of graphene with 
nearest neighbor hopping and is given by 
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where —t < is the hopping integral and ji is the chemical 
potential which can be conveniently controlled by a gate volt- 
age, i^def is the effective Hamiltonian for sp'^ orbitals at sites 
of type 1 
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FIG. 2. (Color online) Plot of the density of states p{E) as a function 
of E — II. By increasing vacancy concentration, /ihaif-fiUed be shifted to 
the left of the Dirac neutrality point as the acceptor band to the left of 
the Dirac point becomes stronger and accommodates more electrons. 



where Esp2 < is the energy level of the sp^ orbital and the 
notation d[^ is adopted for the creation operator at this site 
to emphasize its localized nature in analogy with the situa- 
tion where the localized orbital arises from a transition metal 
ad-atom. Since in this paper we are interested in the physics 
of RKKY interaction, we assume that the local moments are 
formed, and concentrate on their exchange interaction medi- 
ated by the electronic states of the defective host. Therefore 
we do not include the on-site Coulomb repulsion between the 
electrons in this localized orbital. Finally i^hyb is the hy- 
bridization which is given for each defect by: 






(4) 



Here a^cr for i = 2, 3 is an annihilation operator of a pz elec- 
tron at site i with spin a =t,i and V is the amplitude of 
hybridization. The missing carbon atom can belong to both 
sub-lattices. The Jahn-Teller distortion preferring atom 1 over 
the other two can take place in any direction. In this paper 
we take all three possibilities with equal probability. In Eq. |4l 
we set F = — 0.2t for the amplitude of hybridization between 
Pz and sp^ orbitals 111911 and £^^^2 = — 0.5t, unless otherwise 
specified. 

To have a feeling for what happens to the density of states 
(DOS) of the vacancy alloyed graphene, in Fig. [2] the DOS 
of clean honeycomb lattice has been compared to highly va- 
cant graphene. By increasing the concentration of defects, two 
peaks to the left and right of the Dirac point develop which can 
be considered as the analogue of bonding and anti-bonding 
states formed due to the Hybridization term. The chemical 
potential can be tuned with a gate voltage. In the absence 
of a gate voltage and for neutral vacancies we have the half- 
filled situation corresponding to one 2pz electron per existing 
carbon atom. The chemical potential corresponding to half- 



filling deviates to the left of zero when the defects are intro- 
duced to the graphene sample. Therefore the /i = in the 
vacant graphene corresponds to electron doping. 

The general perturbative expression for the RKKY ex- 
change interaction is given by OOiblll . 
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where S is the magnitude of the impurity spin, J is the interac- 
tion between the localized moments and the spin of the itiner- 
ant electrons, i and j are the site index of magnetic impurities 
which are located at position fi and fj, f{uj) is the Fermi- 
Dirac distribution function. Ignoring the constants appearing 
before the integral, the integral appearing in this equation for 
our purposes can be simplified to. 



NUMERICAL RESULTS FOR THE RKKY EXCHANGE 
RKKY interaction in dean graphene 

Let us first discuss the performance of KPM in addressing 
the RKKY interaction. In Ref. IS, the authors used a lattice 
Green's function method to calculate the RKKY interaction 
for pristine graphene beyond the low-energy model of Dirac 
fermions. Their result agrees with other authors flll |23|l. They 
obtain the following analytical results for the clean honey- 
comb lattice: 



^AR — oGJ 



1 



iQ-R 



20r] 



7° 



-CJ' 



[R/af 
cos[(5 • R] 



{R/af 



(7) 
(8) 



Udujf{w)Y^ 



nm 



{En 



■iS){En 



■iS)' 



(5) 



where F^^ = ilj^{fi)ilJn{fj)tlJ^{fj)tlJm{ri), with i/jniu) de- 
noting the eigenvector corresponding to the eigenvalue En of 
the host Hamiltonian. The lattice constant a and h are set to 
unity in all calculations. We assume that temperature is zero 
(T = 0) so that the Lehman representation of the above inte- 
gral in terms of appropriate spectral functions becomes. 
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where C = 9A^fi^/(2567rt) is a positive quantity, a is the lat- 
tice constant, K^ = (±27r/(3\/3a), 27r/(3a)) are the Dirac 
points in the momentum space, R = fi — fj, and angle Or 
is the polar angles between R and Q = i^+ — K~. These 
equations indicating ferromagnetic (antiferromagnetic) corre- 
lation on same (different) sub-lattice were numerically con- 
firmed with the KPM method IlZI, [sl]. We have also checked 
that the KPM results obtained in the clean limit (see Fig. [3]) 
agree with the numerically exact method of Green's functions 
for the pristine graphene DSl lZD. The strength of the KPM 
method is that it can handle clean and disordered systems on 
the same footing. For the disordered case only an averaging 
over realizations of disorder will be required which is what we 
explain in detail in the following subsection. 



where F{e^e^) = Rc[Aji {e)Aij (e^)] and the real-space "spec- 
tral function" is given by Aij (e) = {i\S{e — H)\j). The kernel 
polynomial method (KPM) can be conveniently used to eval- 
uate various spectral functions including A^j (e) ||7l|9l|22]. In 
this method, A^j (e) can be expressed as a series in Cheby- 
shev (or any other complete set of orthogonal) polynomials 
and the expansion coefficients - the so called moments - are 
evaluated by repeated operation with powers of the appropri- 
ate re- scaled Hamiltonian in such a way that a T^ (i^) is gener- 
ated via their recursive relations. In the appendix B we briefly 
summarize the KPM method. The defects considered here 
are on the scale of few percent and their positions in the lat- 
tice as well as the preferred direction due to Jahn-Teller dis- 
tortion is drawn from a uniform random distribution. Typi- 
cal averages are obtained by geometric averages of the form 
^RKKY = exp(((l/2)ln(jRKKY)^)ave). The gcomctric av- 
erages of the above type are known to better represent the 
propagating nature of waves, and when they become zero in- 
dicate the transition to insulating state iAw- To ensure the 
same holds for our particular form of disorder, in appendix A 
we have compared the statistics of geometric and arithmetic 
averages. 
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FIG. 3. (Color online) Plot of the Jrkky as a function as R in the 
armchair direction. By increasing distance R between two magnetic 
impurities an atomic- scale sign alternation of Jrkky agrees with the 
analytic formula |5]. For a 1% concentration of vacancies the oscil- 
lations due to interference between the two valleys are washed out. 
The change in the exponent a depends on the chemical potential. 



RKKY interactions in heavily vacant graphene 

We create evenly distributed defects on the graphene in or- 
der to measure their effect on the mediation of RKKY inter- 
actions. The energies are measures in units of the hopping 
ampHtude, i.e. we have set t = 1. We used a graphene sheet 
with periodic boundary condition and 2.5 x 10^ lattice points 
are used. We generate 1.5 x 10^ different random configura- 
tions for the distribution of vacancies. The number of Cheby- 
shev moments is in the range 3000 — 4500. The hybridization 
parameter is fixed atV = — 0.2t, and Egp2 = —OM. 

In Fig. [3] we have plotted the distance dependence of the 
RKKY interaction for a sample with 1% vacancies. The in- 
set shows DOS and the integrated DOS (electron number) for 
clean and 1% vacant sample. The results are reported for 
the armchair direction The first important feature in pristine 
graphene is the presence of atomic scale oscillations which 
are due to the interference between two valleys. The second 
feature which has become more manifest in the log-log plot 
is that the atomic scale oscillations diminish in long distances 
for both /ihaif- filled and /i = cases. This means that the Bril- 
louin zone scale wave vector connecting the two Dirac cones 
of the pristine graphene does not exist in 1% vacant sample 
anymore. 

Another feature seen in the two plots corresponding to va- 
cant sample is that the slope of the exponent a characteriz- 
ing the overall power-law decay of RKKY interaction stays at 
a = —3 for the half-filled situation, while it is below a = —3 
for the electron doped (/i = 0) case. To further explore this 
observation, we check that when the chemical potential falls 
in the p{E) oc \E\ region (which is naturally inherited from 
the density of states of graphene parent) the exponent remains 
at a = —3 in agreement with other KPM study |[8|]. Indeed a 
linear average DOS in two dimensions corresponds to a linear 
dispersion relation, which by dimensional analysis is expected 
to give rise to R~^ dependence in the particle-hole bubble. 
This explains why in low concentration regime, despite intro- 
ducing 1% vacancies, the exponent a = — 3 in Fig.O With 
this point in mind, we now proceed to study the variations in 
the exponent a in vacant graphene at /i = 0. 



Doping dependence 

Based on the discussion in previous subsection, at /i = 
which corresponds to electron doping in vacant graphene and 
the chemical potential falls in the region of DOS which is 
strongly nonlinear, we expect significant variations in the ex- 
ponent a which can be attributed to hindered motion of elec- 
tron waves inside the vacant graphene medium. For this, in 
Fig. [4] we plot the dependence of RKKY interaction to dis- 
tance in a log-log scale for the armchair direction. To extract 
the exponent a in the clean and vacant systems, first we dis- 
card the length scales below R = 10a, as the power-law is 
meant for long distances. Secondly, slight oscillations surviv- 
ing the presence of disorder are removed in order to remain 
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FIG. 4. (Color online) (a) Plot of the Jrkky as a function as R in 
the electron doped situation (/x = 0). By increasing distance be- 
tween two magnetic impurities in the armchair direction, Jrkky still 
decreases with by a power-law. (b) Plot of the a as a function of 
vacancy concentration in armchair direction. By increasing vacancy 
concentration, a will decrease. For the zigzag direction the trend is 
similar, with slightly different exponents at a given concentration. 



focused on the power-law part of the dependence. We will 
return to the question of oscillations in the sequel. 

The clear change in the exponent a indicates that the prop- 
agation of electrons in heavily vacant graphene is harder than 
the clean sample which is conceivable as the disorder is ex- 
pected to make the propagation of both single-particle and 
particle-hole excitations more difficult. Based on dimensional 
arguments, the change in a could be effectively ascribed to 
a change in the average dispersion relation of the electronic 
states characterizing the vacant sample. Along the zigzag di- 
rection a similar conclusion holds, with a minor difference 
that for a given concentration, the value of the exponent a 
for the zigzag direction slightly differs from the armchair di- 
rection. The qualitative distinction between the armchair and 
zigzag directions of the pristine graphene does not prevail to 
the realm of high concentration of vacancies and the qualita- 
tive behavior along the two directions are not expected to be 
much different when a few percent vacancies are introduced. 

At the /i = situation we also examine the dependence of 
a on the energy level Esp2 of the vacancy. As can be seen in 
Fig. [5] this exponent is not sensitive to the precise value of the 
impurity level as long as it is not deep enough to allow for the 
formation of localized bound states which totally changes the 
nature of propagation of electronic waves. 



The fate of Fermi-surface oscillations 

When the pristine graphene is doped, in addition to atomic 
scale oscillations in the RKKY interaction which merely arise 
from the wave vector Q connecting the two Dirac cone val- 




FIG. 5. (Color online) Plot of the Jrkky as a function as distance for 
(a) armchair and (b) zigzag directions at /i = for different values 
of Egp2 . As can be seen, the exponent a is not much sensitive to the 
precise value of Egp2 . 



leys, there appears another longer range oscillations due to 
the presence of extended Fermi surface with characteristic 
wave- vector kp Mi • When the disorder is introduced within 
the Anderson model, both atomic scale oscillations and the 
Friedel oscillations due to Fermi surface survive [8]. Note 
that the Dirac picture remains quire robust against the weak 
disorder ^ • Therefore given the dimensional consistency the 
overall physics of the RKKY interaction is expected to sur- 
vive in weak disorders considered in |i8J. As we show in 
the following, for the model of vacant graphene considered 
here, even the Fermi surface oscillations are washed out. To 
begin with, note that since the functional dependence of the 
oscillatory part of the RKKY interaction on distance R al- 
ways comes through a combination kpR, it could be viewed 
as an oscillatory function of kr instead. With this in mind, 
in Fig. Owe have plotted the /i-dependence of the pristine and 
vacant graphene with 1% and 6% vacancies at a fixed distance 
R = 32a. For the clean graphene first of all a perfect particle- 
hole symmetry is seen in the RKKY interaction which is a 
symmetry of the Hamiltonian. Secondly an oscillatory depen- 
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FIG. 6. (Color online) Plot of the Jrkky as a function as /i for arm- 
chair direction (R = 32a). Vacancies tend to wash out Friedel oscil- 
lations. More vacancies do a better job. 



Q^ 



0.4n 

0.3 

0.2 

0.1 

0.0 
-0.1 
-0.2 



clean |li=0 
clean |li=0.1 
concentration=1%, ^i^^,^.^.,,^^ 

concentration=1%, |li=0 
concentration=1%, |u=0.3 




20 



40 
R 



60 



FIG. 7. (Color online) Plot of the i?^ Jrkky as a function as R for 
armchair direction. The clean limit is compared with 1% vacancy 
situation at different doping levels. Adding one percent vacancies 
kills both type of short-ranged and long-ranged oscillations for the 
considered chemical potentials. 



dence on /i is seen which implies oscillatory dependence on 
kp and hence on kpR on dimensional grounds. By adding 
vacancies to the system note that the particle-hole symmetry 
is lost as the presence of vacancies breaks such symmetry. 
Now at 0.01 vacancy concentration the there are no Friedel 
oscillations up to /i values where the pristine graphene would 
complete one cycle of oscillations. Some oscillations with re- 
duced amplitude can be seen below the /i ^ —0.3. When 
the vacancy concentration is increased up to 0.06 no Friedel 
oscillations can be found. 

To compare the real-space profile of the oscillations more 
clearly, in Fig. [7] we have plotted the i^^ Jrkky as a func- 
tion of R along the armchair direction for various values of 
/i at vacancy concentration 0.01. The plots are contrasted to 
the case of pristine graphene. As can be seen in the pristine 
graphene, when the chemical potential is at /i = which cor- 
responds to undoped case, the function i^^ Jrkky only shows 
an atomic scale oscillations at large R (black, filled-square 
plot). When the chemical potential is tuned away from the 
Dirac point in pristine graphene (the red, filled-circle plot) it 
can be clearly seen that in addition to the atomic scale oscilla- 
tions, a Friedel type oscillations due to the presence of a Fermi 
surface with a definite characteristic scale kp appear and are 
superimposed into the atomic scale oscillations in agreement 
with analytical results for clean graphene H. However the 
striking observation is made when one percent vacancies are 
introduced. In such case, not only the atomic scale oscilla- 
tions are diminished, but also no sign of oscillations due to 
"Fermi surface" are seen at chemical potentials correspond- 
ing to half-filling and /i = 0.1. At/i = 0.3 it seems that 
some tiny atomic scale oscillations are present, while longer 
range Fermi-surface oscillations are still absent. The absence 
of Friedel-type oscillations in vacant graphene indicates that 



in presence of such amount of vacancies the picture of a Fermi 
surface with a sharp length scale kp^ does not hold anymore. 
It is reminiscent of disorder driven non-Fermi-liquid behavior 
in Kondo alloys 12411 where the creation of vacancies plays a 
role akin to alloying. Therefore it is likely that in such case the 
Fermi surface is characterized with multitude of length scales 
arising from local Fermi surfaces, the disorder averaging over 
which washes out the part of oscillations corresponding to 
presence of a sharp single Fermi wave-vector scale of clean 
metallic state. It should also be noted that in Fig.[7]in the case 
of pristine graphene the amplitude of long-range oscillations 
in the function R^Jrkky increases which indicates the decay 
power in these type of terms is weaker than the R~^ behav- 
ior (fl]. The same argument for the vacant system at /i = 0.3 
indicates that the decay power of the Jrkky is not as strong 
8iS R~^ behavior. 

To summarize we find that when vacancies on the scale of 
a percent are introduced to graphene, the nature of RKKY in- 
teraction changes as follows: when the doping level is such 
that the chemical potential falls in the energy region where 
DOS is linear, the exponent a in R^ remains at a = —3. 
Otherwise presence of vacancies pushes a to more negative 
values. Regarding the oscillatory nature of the RKKY inter- 
action in pristine graphene, a few percent vacancies wash out 
both atomic scale oscillations and the Friedel oscillations due 
to Fermi surface for a remarkable range of chemical potential 
values. 
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APPENDIX 

A: Statistics of geometric averages 

In this appendix we present the statistics of our geometric 
averaging process. We evaluate the probability distributions 
for both X = (1/2) In ( J|kky) (corresponding to geometrical 
averaging) and Jrkky (corresponding to arithmetic averaging) 
in Fig. [8] panels (a) and (b) respectively. As can be seen the 
distribution in x is a Gaussian, 
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By fitting the Gaussian to the distribution data in panel (a), 
we find the corresponding width a of the distribution which 
as been plotted in Fig.[8](c) as a function of the impurity con- 
centration. It is interesting to compare the behavior in Fig.[8]3 
with the linear behavior reported in Ref . H for a the Anderson 
model disorder. 
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FIG. 8. (Color online). Details of the statistics, (a) The distribution 
of the quantity x = (1/2) In ( J|kky) sampled from the armchair 
direction distance R = 60a. (b) The non-Gaussian distribution of 
Jrkky (corresponding to arithmetic averaging), (c) Plot of a in ge- 
ometric averaging as a function of the percentage of vacancies in 
armchair direction. By increasing the concentration of vacancies the 
variance will increase. All data are evaluated for 3 x 10^ configura- 
tions. 



B: Adjusting the chemical potential in KPM 

In this appendix we derive a simple relation that enables us 
to find out the chemical potential corresponding to a given 
number of electrons in the system. To be self contained, 
we briefly review the kernel polynomial method (KPM) as 
weU m. 

Consider a quadratic Hamiltonian H whose energy eigen- 
values E are limited to a bandwidth [EminjEmax]- Rescaling 
the Hamiltonian from H{E) to H{s) where H = {H - h)/a 
and e = {E - b)/a where b = {Emax + Emin)/'^ and 
a = {Emax — £^min)/2; onc can expand spectral functions 
such as the density of states (DOS) in a complete set of e.g. 
Chebyshev polynomials defined in the range e e [—1,1] as, 
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where Tm{£) = cos(marccos(6:)) are the mth Chebyshev 
polynomials, firn are Chebyshev moments and gm are ap- 
propriate attenuation factors to minimize the Gibbs oscilla- 
tions. The moments are traces of polynomials of the Hamil- 
tonian and can be statistically calculated as a trace, /i^ = 
1/r ^r^i {(l)r\Tm{H)\(j)r), whcrc (j)r are a random single par- 
ticle states and M is the number of random realizations used 
in numerical calculations and Nc is a large number where the 
expansion is cut off. To obtain the effect of Tm{H) on a 
given ket, the recurrence relation of Chebyshev polynomials, 
Tm{H) = 2HTm-i{H) - Tm-2{H) with initial conditions 
Ti {H) = H and To{H) = 1 are used. At the end the DOS in 
the original bandwidth scale is obtained as p{E) = p{e)/a. 

The spectral functions corresponding to space non-diagonal 
propagation can also be similarly calculated between any two 
points i and j in the lattice. 
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where the non-local moment /i^ is given by the matrix ele- 
ment (z|T^(i^)|j). 

When we add vacancies in a graphene sheet, the particle- 
hole symmetry is lost, and the half-filling will not correspond 
to /i = anymore. Therefore within the KPM we need to 
find a formula to give a relation between the particle number 
density n and the chemical potential /i. The filling factor is 
given by. 



i(/i) 



/:" 



6{e — fi)p{e)de 



(12) 



where 6{e — ji) is Heaviside function and ji is the (scaled) 
Fermi level. Expanding the integral in Chebyshev polynomi- 
als and using the trigonometric representation of Chebyshev 
polynomials one can easily perform the integration to get, 

n= arccos(— /i) — 2 > sm(marccos/i) (13) 



TT 



^-^ mix 

m—l 



For a given n one needs to adjust /i to satisfy this equation. 
Then the physical /i is obtained by undoing the reseating. For 
example it can be checked that in clean graphene, the solution 
of the above equation with n(/i) = 0.5 will be /i = 0. 
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